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National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 

To study the nonlinear physics of incompressible turbulent flow, the 
unaveraged Navi er-Stokes equations are solved numerically. Initial three- 
dimensional cosine velocity fluctuations and periodic boundary conditions 
are used. No mean gradients are present. The three components of the mean- 
square velocity fl uctuations are equal for the initial conditions chosen. 

The resulting solution shows characteristics of turbulence, such as the non- 
linear excitation of small-scale fluctuations. For the higher Reynolds num- 
bers the initially nonrandom flow develops into an apparently random tur- 
bul ence. 
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INTRODUCTION 


Although much turbulence research has been done in recent years, there 
does not appear to be a general deductive theory of high Reynolds number 
turbulence. Most of the analytical theories depend on a closure assumption 
for the hierarchy of averaged equations. This immediately calls into ques- 
tion the appropriateness of referring to them as deductive theories, except 
at low Reynolds numbers where nonlinear effects are small. ^ 

One way in which a closure assumption can be avoided is by closure by 
specification of initial conditions.^ That method can successfully pre- 
dict turbulent decay. However, in order to use it the initial conditions 
must be fully turbulent, and a large amount of initial data is required to 
satisfactorily specify the initial turbulence. Moreover, the method appears 
limited to decaying turbulence and does not seem capable of extension to 
turbulence maintained by mean gradients, where the effect of initial condi- 
tions eventually becomes negligible. 

In view of the foregoing comments it seems desirable to consider a 
numerical solution of the Navier-Stokes equations which displays features of 
turbulence. It might be pointed out that it is more appropriate to refer to 
a numerical solution of the unaveraged equations as deductive than it is to 
so refer to most of the analytical theories, where closure assumptions are 
required. Several numerical solutions of the uriaveraged equations have 
appeared which use a spectrum of random initial conditions (e.g., Refs. 3 
and 4). These studies appear to demonstrate the feasibility of carrying out 
turbulent solutions with present-day computing equipment and certainly re- 
present major advances. 

Because of the difficulty of specifying realistic turbulent initial 
conditions (experimentally or analytically), it may be more appropriate to 



initially specify a simple regular fluctuation (as actually occurs down- 
stream of a grid). This should be better for studying the development of 
smallscale fluctuations than would a spectrum of initial fluctuations, since 

for the latter, small-scale fluctuations are already present in the initial 
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flow. Taylor and Green and others ’ have used a perturbation series 
to calculate the nonlinear development of higher harmonics from lower ones. 
However, those studies were limited to short times or low Reynolds numbers. 

O 

Orszag and.Fateman have recently used Taylor and Green's initial condi- 
tions and obtained a numerical solution for higher Reynolds numbers and 

larger times. The inviscid (infinite Reynolds number) case was investigated 
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in some detail by Betchov and Szewczyk. 

In the present paper we study the nonlinear physics of turbulence 
numerically. The initial conditions are chosen to give a flow with tur- 
bulentlike properties. Unlike the problem of Taylor and Green, all three of 
the directional components of the mean-square velocity fluctuations are 
equal. This is a simpler type of flow than that considered previously, in 
that it is only necessary to follow the evolution of one velocity compo- 
nent. We are assured that if one component has the desired turbulent or 
turbulent! ike characteristics, the overall flow will have them. That is not 
the case for all initial conditions. For instance initial conditions tried 
earlier gave two mean-square velocity components which were equal and 
decayed rapidly like turbulence, but the third decayed much more slowly. 

One of the objectives is to study the nonlinear excitation of small- 
scale fluctuations (starting with only large-scale ones) for a range of 
Reynolds numbers. This will be done partly by calculating the evolution of 
the microscale of the fluctuations. It will be seen that although the ini- 
tial conditions are nonrandom, the flow at higher Reynolds numbers breaks up 
into an apparently random turbulence. 
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EQUATIONS AND METHODS 

The Navier-Stokes and continuity equations for an incompressible fluid 
can be written in dimensionless form as 


and 


where 


3u i 

W~ 


3 < u i“k> _ i£_ + ’S 


ax. 


ax 


3X k 3X k 



= o. 


( 1 ) 

( 2 ) 



V 




t 


V 

= T 

Xr, 


* 

t , and 


* 



Note that the stars on dimensional quantities are omitted for corresponding 
dimensionless quantities. The subscripts can take on the values 1, 2, 
and 3, and a repeated subscript in a term indicates a summation. The quan- 
tity u i is an instantaneous velocity component, x^ is a space 
coordinate, Xq is a characteristic length, t* is the time, p is the 
density, v is the kinematic viscosity, and p* is the instantaneous 
pressure. In order to obtain an explicit equation for the pressure, we take 
the divergence of Eq. (1) and apply the continuity Eq. (2) to get 
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In the remainder of the analysis it will be found convenient to use tne set 
of Eqs. (T) and (3), rather than Eqs. (1) and (2). 

The expression assumed for the initial disturbance is, in dimensionless 

form 


u^ = a.j cos 
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where 


, _ ^2 a * b - ^ h* c - ^ c* 5 - x ?* f - x f % - x f 

a i " v a i ’ D i " v i’ c i " v c i ’ q “ x O q ’ r ~ x 0 ’ s " x 0 * 

The quantities a^, b^, and are initial velocity amplitudes or 

Fourier coefficients of the disturbance, and q , r , and s are initial 
wave number vectors. The quantities q.x, etc. are dot products, so that 
q.x = q -j x i + q 2 x 2 + ^3X3* etc. In order to satisfy the continuity 
condition, Eq. ( 2 ), 
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For the present work we set 

a i - k ( 2 , 1 , 1 ) , b. = k( 1,2,1 ) , c. = k(l,l,2), 

= (-1,1,1), r. = (1,-1, 1), and $. = (1,1,-1), 


( 5 ) 
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where k is a quantity that fixes tne Reynolds number. In addition to sat- 
isfying continuity, Eqs. (6) give 


2 2 2 T 

U 1 - u 2 - u 3 = u 


( 7 ) 


at all times, where the overbars indicate values averaged over space. Thus, 
Eqs. ( 4 ) and (6) give a particularly simple type of flow, in that we need 
specify only one component of the mean-square-velocity. 

In order to carry out the numerical solution of Eqs. ( 1 ) and ( 3 ) sub- 
ject to initial conditions, Eqs. ( 4 ) arid (6), we use a cubical grid with 
faces at x^ = 0 and 2 *. For boundary conditions we set the values of 
u.j equal at corresponding positions on opposite faces (and in the vicin- 
ity of those faces, as required for calculating numerical derivatives at the 
faces). That is, we assume that the boundary conditions are periodic. 


4 



The spacial- arid time-differencing schemes are essentially those used 
by Clark, et al. 4 That is, for the spatial derivatives in Eqs. (1) and 
(3) we use centered fourth-order difference expressions. 10 For time- 
differencing we use a predictor-corrector method with a second-order (leap- 
frog) predictor and a third-order (Adams-Moulton) corrector. 11 The 
Poisson equation for the pressure (Eq. (3)) is solved directly (no itera- 
tion) by a fast Fourier-transform method. This method of solution was found 
to preserve continuity (v.u»0) except near the ends of some of the runs, 
where the solutions began to deteriorate. (Another indication of incipient 
solution deterioration near the ends of some of the runs was tnat Eqs. (7) 
were no longer accurately satified.) 

Some of the results will be extrapolated to zero mesh size in an effort 
to obtain more accuracy. An expression for a corrected value of the mean 

square velocity fluctuation u£ in terms of the uncorrected value 
~ 

u u , which is consistent with fourth-order spacial differencing, is 

u c = u u + B ( t ) Ax4 + C(t)ax 6 + . . ., (8) 

where ax is the distance between grid points. In a plot of u 2 against 
time t it may be desirable to correct the time for a given u 2 rather 
than u for a given time. For that purpose we have, instead of Eq. (8) 

t c = t u + B(u 2 )ax 4 + C("?)ax 6 + . . . (9) 

In choosing between Eqs. (8) and (9) we use the equation that gives the 
smaller correction. Both expressions are accurate for an infinite number of 
terms, but for a finite number of terms tne expressisn that gives the 
smaller correction should be more nearly converged in the sense that addi- 
tional terms would not significantly change the relation between and t. 
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as a 


The unknown values in Eq. (8) or (9) are obtained by calculating u : 

function of t for several different grid-point spacings ax. For instance 

by truncating Eq. (9) after terms of order 6, we calculate t against 

u^ for three values of ax to obtain t , B(u^) and C(u^) for each 
~Y 

u . If the three values of ax differ successively by factors of two, 

l c = W l l ■ im *2 W5 t 3’ (10) 

where tp tg, and tg are calculated values of t for three values 
of ax at a fixed u . If only two ax's are used. 


t c = TF t l " T5 V (11) 

p- 

Similar equations are obtained for u^ at a given time. 

RESULTS AND DISCUSSION 

Figure 1 shows the calculated evolution of velocity fluctuations at two 
fixed points in space for a high Reynolds number (no correction for dis- 
cretization error). In spite of the nonrandom initial conditions (Eqs. (4) 
and (6)) the velocity fluctuations have the appearance of those for a random 

/-?; 1/2 

turbulence. The dashed curves of Uj/uq are for initial condi- 
tions perturbed ~0.1 percent. The perturbed curves follow the unperturbed 
ones for a short time and then depart sharply. Thus a very small perturba- 
tion of initial conditions causes a large change in the values of u^ 

(except near t = 0). On the other hand, the root-mean-square values of the 
velocities decrease smoothly and are unaffected by the perturbation of ini- 
tial conditions. All of these features are characteristic of turbulence. 

Another characteristic of turbulence is that the correlation between 
velocities at two times should go to zero as the separation of the times 
increases. Figure 2 shows that this occurs for the present high Reynolds 
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number calculations. For true turbulence the correlation should probably 
decrease smoothly with time. Figure 2(b) shows that this is nearly the case 
for the larger t^ (the time that is fixed as the other time increases). 

At early times there is probably some nonrandom structure in the flow caused 
by the nonrandom initial conditions (Fig. 2(a)). 

As a further indication that the high Reynolds number flow breaks up 

into turbulence, we calculate the evolution of the cross correlation 

~~2 ~2 

u^. Although u-^ = u 2 = u^ at all times, the initial UjU 2 given by 

Eqs. (4) and (6) is not zero. However Fig. 3 shows that because of the ap- 
parent randomization of the flow u^u 2 goes to zero as time increases. 

The fluctuations in the curve at early times (as also in the curve of Fig. 
2(a)), are probably caused by nonrandom structure in the flow at early 
times. Figures 1 to 3 show that at later times we appear to get a reason- 
able approximation to isotropic turbulence, although the initial conditions 
are nonrandom. Calculated values of velocity-derivative skewness factor for 
a high Reynolds number (Fig. 4) also appear to be comparable to those for 
isotropic turbulence. The present calculations differ from earlier ones, 
where turbulence was obtained by using random initial conditions. Our 
calculations at lower Reynolds numbers give results that are less turbulent- 
like. Thus the fluctuations develop apparent randomness only at the higher 
Reynolds numbers. 

A question remains as to the source of the randomness observed in 
Fig. 1. Until recently it was generally assumed that randomness in a tur- 
bulent flow is due to randomness in the initial conditions, to random exter- 
nal fluctuations, and/or to the presence of so many eddies or harmonic com- 
ponents (or of so many degrees of freedom) that the identity of the individ- 
ual eddies is lost. * In the present results the first two of these 
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have already been eliminated, thus leaving only the proliferation of eddies 
or harmonic components as a source of apparent randomness. ^ However 
Lorenz and others have shown that a system of nonlinear equations 
with only a few degrees of freedom can develop an apparently random behavior 
in time as a result of loss of stability of the solutions. This happens 
because of the occurrence of so-called strange attractors, which are regions 
in the phase space of the system to which solutions are attracted. 

Randomness arises in those regions from a haphazard movement of the phase 
point among the neighborhoods of various critical points in the phase space 
(points of unstable equiliorium where the solution is steady state). 

Strange attractors can occur even with a small number of degrees of freedom 
(two or three). 

The system of spacial ly-differenced equations used here (disregard time 
differencing), which is obtained from the Navi er-Stokes equations, consists 
of ordinary differential equations of essentially the same form as those in 
the Lorenz system. 14,15 j n particular, both systems have nonlinear alge- 
braic terms that contain products of two dependent variables. Since a 
strange attractor arising from those nonlinear terms is known to produce 
randomization in the Lorenz system, it appears that such attractors should 
also produce randomization in the differenced Navier-Stokes system. For a 
sufficiently fine grid there should be no significant difference between 
results from the undifferenced and differenced versions of the Navier-Stokes 
equations, so that it seems reasonable to assume that strange attractors 
play a role in Navier-Stokes turbulence. The same conclusion is reached if 
the spectral rather than the differenced form of the Navier-Stokes equations 
is used, since those equations also have the same form as the Lorenz 
equations. 
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In an actual turbulent flow, as in the model used here, there is a very 
large number of eddies or degrees of freedom, so that although strange 
attractors probably produce randomization, additional apparent randomization 
will almost certainly be produced by prol iteration of eddies or harmonic 
components. It appears that no conclusions can be drawn as to the relative 
importance of the two processes. The sensitivity of the flow to initial 
conditions and the suddenness with which the effect of perturbation of ini- 
tial conditions becomes manifest (Fig. 1) are results often associated with 

1 01 K 

strange attractors. On the other hand those results could be caused 

in large measure by the nonlinear production of harmonic components, since, 

particularly at higher Reynolds numbers, the latter tends to be explosive 

£ 

(each harmonic component interacts with every other one). 

It might be thought that larger Reynolds numbers would favor randomiza- 
tion by proliferation of harmonic components, since that process is stronger 
at high than at low Reynolds numbers. However the larger number of degrees 
of freedom at high Reynolds numbers will correspond to a greater number of 
unstable equilibrium points in the phase space of the system (where the 
solution is steady state). There will then be more opportunities for 
randomization in strange attractors. Thus at high Reynolds numbers there 
are arguments for randomization both by proliferation of harmonic components 
and in strange attractors, and both will probably occur. 

The presence of strange attractors may be fortunate from a numerical 
standpoint, in that it should enable turbulent solutions that are qualita- 
tively correct (at least insofar as they appear random in time) to be ob- 
tained with a relatively coarse grid. A calculation which used a grid with 
only 3^ points gave results which appeared to be about as random as those 
in Fig. 1, where 32^ points were used. The use of a fine three- 
dimensional grid, of course, requires a large amount of computer time. 
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In the results given so far, no correction for the discretization error 
due to the finite numerical mesh size was applied. The primary purpose of 
the present work is to stuay the physics of turbulence rather than to obtain 
highly accurate results (possibly unattainable at very high Reynolds num- 
bers). For low Reyno las numbers Fig. 5 shows that surprisingly good results 
can be obtained even with coarse grids. At higher Reynolds numbers the 
results, although less accurate, should still be qualitatively correct. 

Their accuracy can be improved by extrapolating to zero mesh size by apply- 
ing fourth-order discretization corrections as explained in the last sec- 
tion. This is done in lieu of subgrid modeling, which has been used in 
previous work.^ The method is related to subgrid modeling in that it 
assumes that the subgrid eddies are somewhat similar to the calculated 
eddies, but does not require the introduction of a subgrid eddy viscosity. 

Iri all of the remaining results the fourth-order discretization corrections 
are applied. However the corrections are negligibly small except at the 
highest Reynolds number. (An alternative (possibly more accurate) method of 
extrapolating to zero mesh size would be to perform the extrapolation at 
each time step (or for a group of steps) as the numerical calculation pro- 
ceeds, rather than after all the calculations for the various inesh sizes are 
completed. However that procedure would require a special program and was 
not used here. ) 

Figure 6 shows the calculated evolution of mean-square velocity fluc- 
tuations for a series of initial Reynolds numbers. As the Reynolds number 

~2 

increases (v and Xq held constant), the rate of decay of u in- 
creases sharply, as in experimental turbulent flows. This can be attri- 
buted to the nonlinear excitation of small-scale turbu lentl ike fluctuations 
at the higher Reynolds numbers. The high shear stresses between the small 
eddies cause a rapid decay. 
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The development of the small-scale eddies is seen more clearly in 
Fig. 7, wltere microscale x calculated from 

*2 

x 2 = _ IQvu * 
du /dt 

(Ref. 16) is plotted against dimensionless time. As the Reynolds number 
increases, the small-scale structure becomes finer. The microscale de- 
creases until the fluctuation level (inertial effect) is low enough that 
viscous forces prevent further decrease. After x decreases to a minimum, 
it begins to grow. This increase of x with time is due to the selective 
annihilation of eddies by viscosity, the small eddies being the first to 
decay. Thus at large times only the big eddies remain. It is this period 

of increasing x that is generally observed experimentally in grid- 

1 ft 

generated turbulence. 

The early period, in which x decreases with time, is of interest as 

illustrative of inter-wavenumber energy transfer. It is often pointed 

out 1 ’ D that turbulent energy is characteristically transferred from big 

eddies to small ones. This is associated with the nonlinear terms in the 

Navi er-Stokes equations and can be tnought of as due to a breakup of big 

eddies into smaller ones, or as a stretching of vortex filaments to smaller 

diameters. In spite of this energy transfer to smaller eddies, the experi- 

1 ft 

mental results generally show a growth of scale. This is because those 
results are usually for the later period shown in Fig. 7 where, although 
energy is transferred to smaller eddies, the annihilation of small eddies by 
viscous action eventually wins out. The early period shown in Fig. 7 is of 
particular interest, in that the nonlinear transfer effects are truly domi- 
nant there; a sharp decrease in scale actually occurs as energy is trans- 
ferred to smaller eddies. 
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Mean-square vorticity or diss i p t at i on fluctuations are plotted in 
Fig. 8, as calculated from Eg. (6.2.9) in Ref. 16. Although the curve for 
zero Reynolds number, where nonlinear effects are absent, decreases mono- 
tonically to zero, the curves for higher Reynolds numbers increase sharply 
for a while and then decrease. Thus the nonlinear terms in the Navier- 
Stokes equations are very effective vorticity generators and greatly enhance 
the dissipation at small ana moderate times. For large times they appear to 
have the opposite effect, evidently because the turoulence itself decays 
rapidly to zero. 

Figure 9 shows mean-square pressure fluctuations plotted against dimen- 
sionless time. The enhancement of the pressure fluctuations, although not 
as great as that of the vorticity and dissipation again appears to be due to 
nonlinear effects. In this case the nonlinear terms on the right side of 
the Poisson equation for the pressure proDably have the greater effect. 

Three measures of the relative importance of inertial (nonlinear) and 
viscous effects are shown for a moderately high initial Reynolds number in 
Fig. 10. The microscale Reynolds number R the ratio of nonlinear 
velocity term to viscous term, and the ratio of pressure to viscous term in 
Eq. (1) are plotted against dimensionless time. The terms are root-mean- 
square values. All of those measures show a variation from a rather iner- 
tial to a weak fluctuating flow. For instance R varies from about 90 

A 

to 0.7. This is a much greater variation than has been obtained experi- 
mentally for a single run. The curves for the term ratios lie somewhat 
below that for R^. They indicate that except at early times the non- 
linear inertial effects associated with velocity and with pressure do not 
differ greatly. 
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CONCLUDING REMARKS 


The results show that, at least at higher Reynolds numbers, an 
apparently random turbulence can develop from nonrandom initial conditions. 
This is indicated by the appearance of the instantaneous velocity fluctua- 
tions and by the sensitivity of those fluctuations (and the insensitivity of 
average values) to small perturbations in the instantaneous initial condi- 
tions. Moreover, the two-time velocity correlation becomes small as the 
time between the occurrence of the two velocities increases. In addition 
the correlation between two components of the velocity becomes small as the 
time increases, as a result of the randomization. This correlation is not 
small initially, even though the three components of the mean square veloc- 
ity fluctuation are equal at early as well as at late times for the initial 
conditions chosen. Also the calculated velocity-derivat i ve skewness factors 
appear to be comparable to those for isotropic turbulence. Thus except in 
the initial period we evidently obtain, at the higher Reynolds numbers, a 
reasonably good approximation to isotropic turbulence. 

Tne source of the observed randomness may lie in the presence of 

1 2 

strange attractors in the phase space of the system, as well as in the 
occurrence of a very large number of eddies or harmonic components (large 
number of degrees of freedom). However it appears that no conclusions can 
be drawn as to the relative importance of the two processes. (A strange 
attractor is a region in the phase space of the system to which solutions 
are attracted and in which randomness arises from a haphazard movement of 
the phase point among the neigborhoods of critical points (points of un- 
stable equilibrium where the solution is steady state). It can occur even 
with a small number of degrees of freedom.) 
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At early times the nonlinear transfer of energy from big eddies to 
small ones is almost completely dominant and causes a sharp decrease in the 
size of the microscale. Tnis has not been generally observed experimentally 
or analytically because the period usually studied is for later times, where 
the anihilation of small eddies by viscous action causes the scale to grow, 
even though energy is being transferred to smaller eddies. This later per- 
iod of scale growth is also observed in the present results. 

The nonlinear terms in the equations of motion are very effective vor- 
ticity generators arid increase the dissipation and the rate of decay. The 
increased rate of decay is related to the nonlinear transfer of energy to 
smaller eddies; the small eddies decay faster than the big ones because of 
the higher shear stresses between the small eddies. Calculation of terms in 
the equations of motion shows, as might be expected, that the flow is domi- 
nated by nonlinear inertial effects at early times and by viscous effects at 
later times. 

The accuracy of the numerical results can be improved somewhat by 
extrapolating to zero mesh size. This is done by applying a fourth-order 
correction for discretization which is consistent with the fourth-order 
spacial differencing used in the calculations. This procedure is used in 
lieu of subgrid modeling and does not require the introduction of a subgrid 
eddy viscosity. 
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Fig. 5 Effect of numerical mesh size on evolution of u^ 
at low and moderate Reynolds numbers. 



Fig. 6 Calculated evolution of mean-square velocity 
fluctuations (normalized by initial value) for various 
initial Reynolds numbers. 
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Fig. 7 Calculated evolution of microscale of velocity 
fluctuations (normalized by initial value) for various 
initial Reynolds numbers. 



Fig. 8 Calculated development of mean-square 
vorticity fluctuations u? or dissipation e 
(normalized by initial value) for various inital 
Reynolds numbers. 
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Fig. 9 Calculated evolution of mean -square pressure 
fluctuation (normalized by initial value) for various 
initial Reynolds numbers. 
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Fig. 10 Three measures of relative importance of iner- 
tial and viscous effects versus dimensionless time. 
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